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Operational Tau method for nonlinear 
multi-order FDEs 


P. Mokhtary 


Abstract 


This paper presents an operational formulation of the Tau method based 
upon orthogonal polynomials by using a reduced set of matrix operations for 
the numerical solution of nonlinear multi-order fractional differential equa- 
tions(FDEs). The main characteristic behind the approach using this tech- 
nique is that it reduces such problems to those of solving a system of nonlinear 
algebraic equations. Some numerical examples are provided to demonstrate 
the validity and applicability of the method. 


Keywords: Fractional differential equations(FDEs); Caputo derivative; Op- 
erational Tau method. 


1 Introduction 


The mathematical modelling and simulation of systems and processes based 
upon the description of their properties in terms of fractional derivatives, 
naturally leads to differential equations of fractional order and to the necessity 
to solve such equations. However, effective general methods for solving them 
can not be found even in the most useful works on fractional derivatives and 
integrals. 

There are several approaches to the generalization of the notation of dif- 
ferentiation to fractional orders e.g., Riemann-Liouville, Grunwald-Letnikov 
and Caputo. We focus on one particular form so-called Caputo derivative. 

Recently, linear FDEs based upon the fractional derivatives(such as 
Riemann-Liouville and Caputo schemes) with general variable coefficients 
have been solved by adapting various analytical and numerical methods[2, 5, 
7, 20]. Nowadays, applications have included some classes of nonlinear FDEs, 
and this motivates us to consider their effective numerical methods for solu- 
tion of these type of equations. Among the most recent works concerned with 
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nonlinear initial value problems of fractional order, we can consider papers 
[4, 8, 12, 13, 18, 19, 21, 22, 31, 32]. 

Spectral methods have been studied intensively in the last two decades 
because of their good approximation properties. The formulation of spectral 
methods was first presented in the monograph of Gottlieb and Orszag [11]. 
The text book of Canuto, et al [3] focuses on practical and theoretical aspects 
of global spectral methods. 


Global spectral methods use a representation of function u(t) throughout 
the domain via a truncated series expansion with suitable basis functions. 
This series is then substituted into functional equation and upon the mini- 
mization of the residual function the unknown coefficients are computed. 


Spectral methods can be broadly classified into three categories, Pseu- 
dospectral or Collocation, Galerkin and Tau methods. The Tau method, 
through which the spectral methods, as shown in [6, 23-29] has found exten- 
sive application for the numerical solution of many operator equations in the 
recent years. 


The Tau method, firstly introduced by Lanczos[15-17], involves the pro- 
jection of the residual function on the span of some appropriate set of basis 
functions, typically arising as the eigenfunctions of a singular Sturm-Liouville 
problem. The auxiliary conditions imposed as constraints on the expansion 
coefficients. It is well known that eigenfunctions of certain singular Sturm- 
Liouville problems allow the approximation of functions belong to the space 
C@™/a,b] whose truncation error approaches zero faster than any negative 
power of the number of basis functions used in approximation, as that num- 
ber(order of truncation N)tends to oo. This phenomenon is usually referred 
to as ”Spectral accuracy” (Gottlieb and Orszag [11]). A convergence analy- 
sis and error bounds for the Tau method was considered by Ortiz and Pham 
in the papers [24, 25]. The recursive form of the Tau method, formulated 
by Ortiz in [26] was extended to the case of systems of ordinary differen- 
tial equations in [6]. The basic philosophy of the method was extended to 
the numerical solutions of the linear and nonlinear initial value, boundary 
value, and mixed problems for ordinary differential equations [23, 25, 27], to 
the eigenvalue problems [27, 28], to the ” Stiff’ problems [23], to the partial 
differential equation[29], among others. 


The main objective of the present paper is to provide Ortiz and Samara|27| 
operational approach to the Tau method for the numerical solution of non- 
linear FDEs of the general form 


Lp(u(t)) = f(t), (1) 
on t € A = [0,1] with initial conditions 
w(0)=d;, i=0,1,...,¥—-1, (2) 


where 
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Na Up 
Lp(u(t)) = >— pr(t) [] (Dee ult)", Ore Qt, Nav Yreslr €NLJ{O}, (3) 
r=0 k=0 


where N, Q* are the collections of the all natural and positive rational num- 
bers, respectively. d; are constants and v = ognax, { feild tae The symbol 
r<Na 


[q] is the smallest integer greater than or equal to qg. u(t) is unknown func- 
tion, p-(t) and f(t) are algebraic polynomials or their suitable polynomial 
approximations. Finally, the fractional derivative is considered in the Ca- 
puto sense that is given by 


t 
1 
pout) = p— rial Pray TaD (pdr, EAL (A 
ul) = soggy f to ulnD(r)dr, ted. (A) 
0 


The properties of Caputo derivative can be found in [30]. 


In this paper we proceed as follows: In the next section, the spectral 
Tau method for nonlinear FDEs is described. We reduce the problem to a 
set of nonlinear algebraic equations using some useful operational matrices. 
Numerical experiments are carried out in Section 3, to illustrate the efficiency 
of the proposed method. 


2 Numerical approach 


Consider the operational Tau solution for nonlinear FDE (2- 3) as a polyno- 
mial of degree N 


un (t) = ude) =un J=unJXt, (5) 
i=0 
where uy = [uo,u1,-..,un,0,...]. J is non-singular lower triangular coeffi- 


cient matrix given by the shifted Jacobi polynomials in A, where {J;* Be) 24 
J = JX, with a standard basis vector X; = [1,t,t?,...]"[3]. The effect of 


u) (t), t?un(t) and (uy (t))? on the coefficients vector of polynomial (5) are 


u(t) = unIn*X, tun (t) = unde X:, (un(t))? = un JEP (uy, J) Xt, 


where matrices 7 and yz have the following simple structures ({27]) 
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000. a ee Ome 
10 0 001 9: 
n=|02 0 me a os ee |e 
90 3 0 90 


0 9 


and E(uw, J) is an infinite upper triangular Toeplitz matrix with the following 


structure 
[ edo ads tad 
0 unJo und. ... 
B(uv,J)=| 0 0 undo... |> 


where J; is the i-th column of the matrix J. Details for formulation of the 
matrix E(uy,J) can be found in [10]. 

Now, we intend to explain details of the structure of the operational ap- 
proach to the Tau method with Jacobi polynomial bases for the numerical 
solution of the nonlinear multi-order FDEs. Firstly, in the Lemma 2.1, we 
will show that the effect of Caputo fractional derivative DE (un (t)), will 
be represented as the product of a matrix and a vector. Secondly, in the 
Lemma 2.2, we will prove that the product of polynomials can be written as 
the product of a matrix and a vector. Finally, in the Theorem 2.3, we will 
give the matrix representation of Lp(un(t)) by using the Lemmas 2.1 and 2.2. 


Lemma 2.1 Let Tee) be the shifted Jacobi polynomials with respect to 
the weight function y%*(t) = (2 — 2t)*(2t)? on A. Assume that the ap- 
proxzimated solution un (t) and the fractional derivative Des are given by the 
relations (5) and (4) respectively, then 


De* (un (t)) = un Yo,,J, 


where ; ; 
A - A 
O([Ork | )Em,.p,,0 Of Ore )Epa,.p1.N 
Vor, = | O(fre141)€re,g141.0 - Ore 1+ Epeg ean 2 Ye 


O(N )En,o i O(N)En.n 


with (6) = reg and 
1: k- a,B k > Orr 
| I, JOP te, 
TRO ag FOE GH Oty 
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Proof. See [9]. 


Lemma 2.2 (a) For two given polynomials h(t) = >> hivi(t) = HVXt 
i=0 


and s(t) = > 8;u;(t) = SVX; with H = (ho, hi, he,.--],S = [S0, $1, $2,---], 
i=0 


we have 
s(t)h(t) = SVE(A, V) Xt 


ee 


(b) For given polynomials h;(t) = }> ajTij(t) = awTi Xt, 1 = 0,1,. 
j=0 


where T; are nonsingular coefficients matrices given by {Ti;}72o = TiXz, we 


have 
l 
[[ til) = aw Do [[ Blew, TX. (7) 


Proof. For proof of part(a) see [10]. By using part (a) and mathematical 
induction we can prove part (b). 


Theorem 2.3 (Matrix representation for nonlinear part) 
Assume that the approximated solution uy(t) and the nonlinear fractional 
operator Lp are given by the relations (5) and (2), respectively, then 


Lp(un(t)) = unl J, 


where 
Na L, 
=I (D2 WE (wy, IVa) II E(un, Fre) pr(qt)) Io}, 
re0) k=d4+1 
One 0 N 
— 1) ’ rk = > _ waca 
we {ey Or~ € Qt —N, Fre = IVE (un, JVx), 


and d is the smallest index that 7a 4 0. 
Proof. From Lemma 2.1 and the third relation in (6) for k 4 {p |y¥rp = 0} we 
have 


L,. L,. L,. 


J] (2e* un)" = TT ed XL) = T] (wT WE (uy, IV) XY). 
k=0 k=0 k=0 


Let d be the smallest index that 7,q 4 0, then from (6) we can write 
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L,. 
Il (un JU, EM! (un, JW;,) Xz) =unJV gE! (un, JW.) 


k=0 
L, 
« [J Eun, JU.E™ (uy, JV,))X, 
k=d+1 
= unll,J, 
where 


Ly. 
IL, = JU jE"! (un, IVa) II (un, F,,)J~‘andF,, = JU, E71 (uy, JU,). 
k=d+1 


By substituting the above relation in (3) and using the second relation in 
(6) we obtain 


Na 
Lp(un(t)) = wy (So Trpr(u)) J, 
r=0 


that is the statement of the Theorem. 


Also for obtaining the matrix form of the initial conditions (3), we intro- 
duce vector d = [do,di,...,d)-1,0,...] where vy = max Ogi On the 
0<r<Na 


other hand we can write 


uf) (0) = unIn Xs. = unIner = undj, 7 =0,1,2,...,7-1, 


|.=0 


where e; = [1,0,0,...]7 and 
B = (b3)j=9 = (In'e1)j=5. (8) 


Now, we are ready to obtain the nonlinear algebraic system of implement- 
ing the operational Tau method to the nonlinear multi-order FDE (2-3). 


Following Theorem 2.3 and the relation (7) we obtain 
(9) 
unB=d, 


where f(t) = fJ with f = [fo, fi, ...]. Because of orthogonality of {JP ) hay 
projecting (9) on the {2° (t)}4_, yields 


unlly = fe, = 0,1,2,...,.N. 
By setting 


Mn = (Bo; Biy ++ +3-Dp—44 Lo, Ihh,. ++, lly], en = (do, di,...,dv-1, fo, fi,---, fn], 
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we obtain uy My = ry. We restrict this system to its first N + 1 columns. 
The square system un My_, =Trn-_ p, gives us unknown vector uy. 


3 Numerical results 


In this section we have considered three test problems. All of these test prob- 
lems have been solved by the operational Tau method based on the Chebyshev 
and Legendre bases. In all cases any non-polynomial functions were replaced 
by a suitable polynomial approximation. All calculations were performed on 
a PC running Mathematica software. To report some information about the 
number of operations, we use function \LeafCount in the Mathematica soft- 
ware, that gives us total number of indivisible subexpressions. All of achieved 
nonlinear algebraic systems were solved by the well known iterative Newton 
method. 


Example 1: /18] Consider the nonlinear FDE with a = 1.5 


DSu(t) + u(t) = f(t), u(0)=0, i=0,1, 


[(6) 5—a 


- 305) ja, 20(4) 
T(6— a) T(5 — a) 


p-% + (4° — 3t4 + 203)?. 
T(4— a) ( en) 


the exact solution is u(t) = t? — 3t4 + 2¢°. 


We apply the proposed operational Tau method to obtain the approx- 
imated solution of the problem. The maximal error, with the Chebyshev 
and Legendre bases have been given in Table 1. A comparison of the Tau 
method with fractional high order method proposed by R. Lin and F. Liu 
in [18] shows that our method produces powerful superiority with respect to 
the proposed method in [18]. 


Table 1: Numerical results of Example 1, using operational Tau method with different 


bases 
Maximal Error 
N Chebyshev Tau Legendre Tau 
5 458 10-1 Toe LOL 
7 1.50 x 10716 7.67 x 1071” 
10 8:75: 10°*" 55210. 


Lin and Liu scheme Max. error is 1.544 x 1075 int =1 
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Example 2: Consider the nonlinear FDE 
1 1 - Rest 
Déu(t)D2u(t) + u(t) = et + eerf(V4)(1— as 
4 


where erf(z) gives the error function and T(a, z) is the incomplete Gamma 
function. 


F u(0) =1, 


the exact solution of this problem is u(t) = e*. We apply the proposed 
operational Tau method to obtain the approximate solution of this example. 
We have reported the obtained numerical results in Table 2 and Fig. 1 with 
the Chebyshev and Legendre bases. In Fig. 1, obtained numerical errors 
are plotted for several values of approximation degree N in the LZ. norm. 
From Table 2 and Fig. 1, we can conclude that desired spectral accuracy 
is obtained for this nonlinear problem and the approximate solutions are in 
high agreement with the exact solution. In addition, by using the function 
\LeafCount in the Mathematica software, for N = 4,8,12 and N = 16, we 
need 213,812,1261 and 2153 operations, respectively, to obtain the opera- 
tional Tau solution with the reported errors in the Table 2 and Fig. 1 based 
on the Chebyshev polynomial bases. 


Table 2: Numerical results of Example 2, using operational Tau method with different 


bases 
Maximal Error 
N Chebyshev Tau Legendre Tau 
4 6.79 x 10-5 bald xel0e® 
8 7.19 x 1071! 5.85 x 1071! 
12 7.31 x 10716 7.06 x 10716 
16 2.85 x 10716 2.81 x 10716 


Example 3: /14] Consider the following equation of fractional order 6 = 0.5: 
Déu(t) = At? (u(t))?, (0< <1), (10) 
with A,B E R(A 40). If64+ 8 <1, this equation has the exact solution 


wty= RECA en a 


If 6 < —20, then the equation (10), has unique solution u(t) € C[a, bj 
given by (11) and if 6 = —(k+0),k EN, then the equation (10) , has unique 
solution u(t) € C™[a, 6]. (See [14, Chapter 3]). 
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Log, o( Error) 
Log, 9( Error) 


Figure 1: An illustration of the rate of convergence for the Tau method with various 


N. We observe the errors of Example 2 using Chebyshev bases (left) and Legendre bases 
(right) 


The numerical results for example 3 with the Chebyshev and Legendre 
bases are presented in Fig. 2 and Table 3. Fig. 2, shows the rate of conver- 
gence for various (. Each part of the figure contains numerical errors for sev- 
eral values of N, which are plotted for a special value of 8 in L., norm. As we 
can see from Table 3., and Fig. 2, the performance of the spectral Tau method 
with the Chebyshev and Legendre bases for 6 € [—1.5,—2.5] almost same, 
but when £ tends to the 6 = —2.5(smooth solution) the rate of convergence 
increases and we have accurate numerical solutions. For 6 = —1.5, —2.5, nu- 
merical results have not been presented, since the exact solution is obtained. 
In addition, by using the function \LeafCount in the Mathematica software, 
for N = 15 we need 917 operations to obtain the operational Tau solution 
based on the Legendre polynomial bases. 


Table 3: The numerical results of Example 3 with different 8, \= 1 and N = 15 


Maximal Error for Chebyshev Tau Maximal Error for Legendre Tau 
x B=-—2 B= —2.2 B=—24 B=-2 B= —2.2 B=—24 
0.2 9.92x10-® 3.38x10-% 5.74x 10-7 1.29x 10> 4.27x 10-5 7.01 x 10-7 
0.4 7.45x10-§ 2.57x10-§ 4.41 x 107-7 1.02x 10-5 3.46 10-8 5.80 x 10-7 
0.6 439x108 1.44x10-° 2.38 x 10-7 7.77x10-° 2.61x10-% 4.31 x 10-7 
0.8 3.47x10-& 1.12x 107° 1.82 x 1077 5.41x10-® 1.70x10-° 2.63 x 10-7 
1 2.00x10-° 5.10107" 9.23 x 10-8 5.20x 10-° 1.67x10-° 2.64 10-7 
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Log,,( Error) 
Log,,( Error) 


Figure 2: An illustration of the rate of convergence for the Tau method with various 


3B. We observe the errors of Example 3 using Chebyshev bases (left) and Legendre bases 
(right) 


4 Conclusion 


In this paper, we presented a numerical scheme for solving nonlinear multi- 
order fractional differential equations. The operational Tau method was em- 
ployed. Also, several test problems were used to show the applicability and 
efficiency of the method. The obtained results indicate that the new approach 
can solve the problem effectively. 
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